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ABSTRACT 

The covariance matrices of power-spectrum (P(fc)) measurements from galaxy sur¬ 
veys are difficult to compute theoretically. The current best practice is to estimate 
covariance matrices by computing a sample covariance of a large number of mock cat¬ 
alogues. The next generation of galaxy surveys will require thousands of large volume 
mocks to determine the covariance matrices to desired accuracy. The errors in the 
inverse covariance matrix are larger and scale with the number of P{k) bins, making 
the problem even more acute. We develop a method of estimating covariance matrices 
using a theoretically justified, few-parameter model, calibrated with mock catalogues. 
Using a set of 600 BOSS DRll mock catalogues, we show that a seven parameter 
model is sufficient to ht the covariance matrix of BOSS DRll P{k) measurements. 
The covariance computed with this method is better than the sample covariance at 
any number of mocks and only ~100 mocks are required for it to fully converge and 
the inverse covariance matrix converges at the same rate. This method should work 
equally well for the next generation of galaxy surveys, although a demand for higher 
accuracy may require adding extra parameters to the htting function. 

Key words: methods: data analysis - galaxies: statistics - cosmological parameters 
- large-scale structure of the Universe 


1 INTRODUCTION 

The covariance matrix and inverse covariance matrix of the 
band averaged power spectrum are crucial for parameter es- 
timatiou from cosmological spectroscopic surveys. Having 
an accurate estimate of the covariance matrix, and there¬ 
fore the inverse covariance matrix, is paramount in order to 
be able to assign reliable uncertainties in estimated parame¬ 
ters (Percival et al. 2014). Most studies achieve this by using 
a large number of mock samples (Cole et al. 2005; Reid et al. 
2010; Manera et al. 2013, 2015; Anderson et al. 2014; Gil- 
Man'n et al. 2015) setup to match the characteristics of the 
particular survey, and then running them through the data 
pipeline to estimate the covariance matrix via 

E {p; - i^j), (1) 

s 

where N is the number of mocks, 

= ( 2 ) 
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and Pi = P(ki). 

The elements of the sample covariance matrix converge 
as O to their true values, while the inverse covari¬ 

ance matrix elements converge as O {N^/N), where Ab is the 
number of P{k) bins (see e.g., Anderson 2003). This inac¬ 
curacy propagates into derived cosmological parameters and 
inflates their errorbars by a factor of O (1 -I- Ab/A) (Hartlap 
et al. 2007; Taylor et al. 2013; Percival et al. 2014; Dodel- 
son & Schneider 2013; Taylor & Joachimi 2014). Percival 
et al. (2014) found that in order for this extra variance to 
be sub-dominant for Baryon Oscillation Spectroscopic Sur¬ 
vey (BOSS; Eisenstein et al. 2011; Dawson et al. 2013) Data 
Release (DR) 9 measurements ~ 600 mock catalogues were 
needed. The next generation of galaxy surveys will have 
more stringent requirements on the precision of the covari¬ 
ance matrix and will require even larger sets of mock cata¬ 
logues to compute the sample variance. Taylor et al. (2013) 
estimated that up to 10® mocks may be required for the joint 
analysis of future galaxy clustering and weak lensing data. 

There are a number of methods with which one can gen¬ 
erate these mock catalogues that tend to reduce the compu¬ 
tational cost, such as the log-normal method (Coles & Jones 
1991), pinpointing orbit-crossing collapsed heirarchical ob- 
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jects (PINOCCHIO; Monaco et al. 2002, 2013), comoving la- 
grangian acceleration simulations (cOLA; Tassev et al. 2013), 
perturbation theory catalogue generator of halo and galaxy 
distributions method (patchy; Kitaura et al. 2014, 2015), 
pertabuation theory halos method (pthaloS; Scoccimarro 
& Sheth 2002; Manera et al. 2013, 2015), or quick particle 
mesh method (qpm; White et al. 2014). Chuang et al. (2015) 
provide detailed descriptions and a comparison of the effec¬ 
tiveness of mocks generated using these techniques, conclud¬ 
ing that the more efficient approximate solvers can be used 
to reach a few per cent accuracy for clustering statistics on 
scales of interest for large-scale structure analyses. 

However, with surveys such as the Dark Energy Sur¬ 
vey (DES; Frieman & Dark Energy Survey Collaboration 
2013), the upcoming Dark Energy Spectroscopic Instrument 
survey (DESI; Schlegel et al. 2011; Levi et al. 2013), Ex¬ 
tended Baryon Oscillation Spectroscopic Survey (eBOSS; 
Schlegel et al. 2009), Large Synoptic Survey Telescope sur¬ 
veys (LSST; LSST Science Collaboration et al. 2009), Euclid 
satellite mission surveys (Laureijs et al. 2011), and numer¬ 
ous others increasing the volumes to be analysed, the mocks 
must also increase in volume. This can still lead to situations 
where the computational costs of generating the necessary 
number of mocks becomes prohibitive, making it desirable 
to have a method of estimating the true covariance matrix 
using fewer mock catalogues. 

The use of mock catalogues can be completely bypassed 
by looking at the intrinsic scatter of P{k) measurements 
within the sample e.g. using the jack-knife or bootstrap 
methods. The jack-knife method is limited by the fact that 
to build up better statistics one must divide the survey data 
into smaller and smaller subvolumes, limiting the maximum 
scale for which the covariance can be reliably measured (Nor- 
berg et al. 2009; Beutler et al. 2011). The bootstrap method 
performs better but is still limited by the number of subvol¬ 
umes that can be created from the data (see Norberg et al. 
2009, for a detailed discussion of the two methods). 

In principle, the covariance of P{k) measurements can 
be computed theoretically. Nonlinear effects in structure 
growth and highly nontrivial survey windows make this kind 
of computation difficult in practice. Despite this, several re¬ 
cent works demonstrated that this approach can be used 
to derive reasonably good approximations to the covariance 
matrix (for recent work see e.g., de Putter et al. 2012, and 
references therein). Similar efforts have been applied to the 
correlation function (inverse Fourier transform of P{k)) co- 
variance matrices (see e.g., Xu et al. 2012). 

Alternative approaches include using a shrinkage esti¬ 
mation (Pope & Szapudi 2008; de la Torre et al. 2013), 
covariance tapering (Paz & Sanchez 2015) and using a 
small number of mocks while ‘resampling’ large-scale Fourier 
modes (Schneider et al. 2011). 

In this paper we propose a new approach to estimating 
P{k) covariance matrices. We start with a brief overview 
of the theory behind P{k) covariance matrices in section 2. 
Then we describe the mock catalogues along with our pro¬ 
cedure for estimating the power spectrum, true covariance 
and inverse covariance matrix, and their associated uncer¬ 
tainties, from those catalogues in section 3. 

In section 4, we choose a theoretically justified func¬ 
tional form with a small number of free parameters to de¬ 
scribe the covariance matrix and use mock catalogues to 


calibrate numerical values of parameters. Unlike previous 
approaches based on theoretical modelling we put signifi¬ 
cantly less stress (and effort) into computing the actual co- 
variance matrix elements; ‘Back of the envelope’ theoretical 
considerations are only used as a rough guide in justifying 
the functional form and the actual numbers come purely 
from the fit to the mock sample covariance matrices. In sec¬ 
tion 5 we show that a simple seven parameter model is good 
enough to describe the covariance matrix of the BOSS DRll 
sample as computed from a sample of 600 mocks. In the 
range of scales relevant for the baryon acoustic oscillation 
(BAO) peak and the redshift-space distortion measurements 
{0 < k < 0.2h Mpc) our fitting function works exceptionally 
well. 

We also show in section 5 that the procedure converges 
much better than the sample covariance with the number of 
mocks used. At any number of mocks, the htted covariance 
matrix is closer to the final result than the corresponding 
sample covariance matrix, and at ~ 100 the fitted covari¬ 
ance matrix is already statistically indistinguishable from 
the sample covariance matrix computed with N = 600. The 
inverse covariance matrix converges at the same rate as the 
covariance matrix. Demand for higher accuracy may require 
introducing additional free parameters into the fitting func¬ 
tion, but there is no reason why this method should not work 
equally well for the future galaxy surveys. These conclusions 
are summarized in section 6, where we also discuss planned 
further work. 


2 THEORETICAL P{K) COVARIANCE 


For a Gaussian field in a large uniform volume the covariance 
matrix of P{k) estimated in bins of width Sk is 


C\j = 


{2nf {P, + n-^f 
V 2nk^6k 


( 3 ) 


where V is the volume of the survey and n is the number 
density of galaxies (Feldman et al. 1994; Tegmark 1997). 
When the number density of galaxies is not constant across 
the survey this changes to 


C.j = 


(271)^ Pf ~ 
V,B{ki) 2nkfSk 


( 4 ) 


with 



[Pln{r) + 1]^ 
Pfn?{r) 


( 5 ) 


If, in addition, the width of the k-bins is comparable to 1/V, 
the finite volume will result in the coupling of neighbour¬ 
ing P{k) measurements and the Kronecker delta function in 
Eq. (4) will turn into 

=y"dVe-4'“--'=7)U (6) 

V 


The observed volume V is usually highly nontrivial which 
makes the effective volume difficult to compute. Nonlinear 
effects in structure growth and galaxy biasing further com¬ 
plicate matters, adding terms proportional to the bin aver¬ 
aged trispectrum, T{ki,kj), to the off-diagonal elements of 
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Figure 1. The average of the power spectra from the 600 CMASS 
NGC PTHALOS mocks. The error bars are the square root of the 
diagonal elements of the sample covariance matrix calculated from 
all 600 mocks. 



the covariance matrix, resulting in 


(27t)^ 

Kff(fei) 


( p? 

\2nkf6k 


+ T{ki, kj 


( 7 ) 


(see Scoccimarro et al. 1999 and Bernardeau et al. 2002 for 
details). 


3 MEASURING P{K) COVARIANCE FROM 
BOSS DRll MOCKS 


While it is possible, in principle, to compute the covariance 
of P{k) measurements from Eq. (7), highly nontrivial sur¬ 
vey windows and the introduction of terms dependent on the 
trispectrum from nonlinear structure growth make this diffi¬ 
cult in practice. Therefore, in order to obtain an estimate of 
the true covariance matrix, we use 600 BOSS DRll pthalos 
mock catalogues (Manera et al. 2013) to compute the sample 
covariance of the spherically averaged P{k). For simplicity 
we only use the mocks for the North Galactic Cap (NGC). 
We use the same estimator, weighting scheme, and shot noise 
subtraction method as the latest official DRll analyses pa¬ 
pers (see e.g., Gil-Marin et al. 2015). We estimate the spheri¬ 
cally averaged P{k) in 23 bins of width Ak = O.OOSh Mpc“^ 
in the wavelength range 0.0 ^ A: 0.184h Mpc“^. We then 
compute the sample covariance matrix using Eqs. (1) and 
( 2 ). 

Figure 1 shows the average P{k) with the errorbars com¬ 
puted by taking a square root of diagonal elements of the 
covariance matrix Figure 2 shows the elements of 

the reduced covariance matrix defined as 


Ci- 


\J CiiCj' 


( 8 ) 


Assuming that the distribution of individual P{k) mea¬ 
surements is close to Gaussian, the measured Cij follow the 
Wishart distribution (Anderson 2003). In the limit of large 
N the Wishart distribution tends to a Gaussian distribution 
with the covariance matrix 

(CijCki) = [r.krji +rarjk), (9) 


Figure 2. The reduced covariance matrix calculated with all 600 
NGC PTHALOS mocks. See the online article for a colour version 
of this plot. 


where (Ji and Vij are the unknown true variance and cross¬ 
correlation coefficients of power spectrum band estimates. 
The variance in both diagonal, 



( 10 ) 


and off-diagonal, 




( 11 ) 


elements of the covariance matrix, estimated using Eq. (1), 
scale with the number of mocks, N. The cross-correlation 
between estimates of different Cij elements is of order of pj 
and can be safely ignored as the measured rij are of order 
of 0.01 or less. 

Of course, the inverse covariance matrix is the quantity 
of interest for parameter estimation. While we are using an 
unbiased estimator to determine our sample covariance, the 
inverse of the sample covariance will not, in general, be unbi¬ 
ased (see Hartlap et al. 2007, for details). In order to obtain 
an unbiased estimate of the true inverse covariance matrix 
it is necessary to apply a correction of the form 


T = 


A-iVb-2^_i 

N-1 


( 12 ) 


where Nh is the number of P{k) bins (Hartlap et al. 2007). 
The variance in the elements of this unbiased inverse covari¬ 
ance matrix can be found as (see e.g., Taylor et al. 2013, 
Percival et al. 2014 and references therein) 


(ATij A4 'm) = ATy -b B{^ik-^ji -b (13) 


where 

"" “(iv-^iVb -i)(N^^)Vb^-Ty’ 

(A - iVb - 2) 

(AT _ ATj, _ l)(Ar - Ab - 4) ■ 

This leads to a variance in the diagonal elements of 

{AT?,) = (A + 2B)tL 


(14) 


(15) 
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and the off-diagonal elements of 

= (4 + (16) 

To summarize, we estimate a 23 by 23 covariance ma¬ 
trix of P(fc) measurements using 600 mocks. The errors on 
the 276 independent elements of the sample covariance ma¬ 
trix are given by Eqs. (10) and (11) with negligible cross¬ 
correlation. In section 4 we will use the sample Cij elements 
and their errorbars estimated in this way to calibrate the 
parameters of the theoretical covariance matrix. Addition¬ 
ally, we estimate the inverse sample covariance matrix using 
Eq. (12), and obtain errors on the independent elements with 
Eqs. (15) and (16). We will use these in section 5, to compare 
how the sample and theoretical inverse covariance matrices 
converge. 



4 CALIBRATING PARAMETERS OF 

THEORETICAL COVARIANCE MATRIX 


Nontrivial survey volume and difficulties inherent in the the¬ 
ory of cosmological structure growth in the nonlinear regime 
make direct computation of the covariance matrix in Eq. (7) 
a highly nontrivial task. Much effort has been put into un¬ 
derstanding the trispectrum of the galaxy field in translin- 
ear and nonlinear limits (see e.g. Scoccimarro & Frieman 
1999; Sefusatti & Scoccimarro 2005). Here we will attempt 
to construct a relatively simple, few-parameter function to 
approximate the true covariance matrix. This fitting func¬ 
tion will, of course, be a very crude approximation to the 
true structure of the trispectrum. However, for the BOSS 
DRll mocks used here, it seems to achieve desirable accu¬ 
racy over a wavelength range relevant to BAO analysis. 

We start by defining two functions, f{k) - to describe 
the behaviour of the diagonal elements, and g(k) - to de¬ 
scribe the behaviour of the off-diagonal elements in the 
correlation matrix. These functions are defined through 
Cii{ki) = Pif^{ki) and rij = g{ki — kj). The first equation 
is the definition of f{ki), while the second equation implies 
that the reduced covariance matrix depends only on the dif¬ 
ference between the centres of bins, an assumption that, in 
general, does not have to hold. 

f{k) is a fractional error in the P{k) measurement and 
since the sample is weighted in such a way as to optimize 
the P{k) measurement at BAO scales we expect it to be a 
smooth function with a minimum around those scales. We 
adopt a three-parameter function 

f{k) = (afc)V^ (17) 

which we justify later in this section. 

The off-diagonal elements are generated by the window 
function effects (I4fr) and the trispectrum (T). For a sim¬ 
ple case of a uniform sample within a cubic volume and no 
additional selection effects the cross-correlation is 


Sij —>■ g{Ak) oc sinc(u;Afe), (18) 

where oj = L/27t, L is the size of the cube, and 
sinc( 2 ;) = sin(a;)/a:. Non-linear gravitational effects will in¬ 
duce some coupling of fc-modes near the diagonal (Meiksin 
& White 1999; Scoccimarro et al. 1999; Sefusatti et al. 2006) 
which we model by a Lorentzian function 
.,2 


g(Afc) oc 


7 


(Afc)2 -I- 7^ 


(19) 


Figure 3. The scaled diagonal elements of the covariance ma¬ 
trix determined from all 600 mocks as a function of k. With the 
chosen scaling this is the fractional uncertainty of the power spec¬ 
trum. Our fitting function (/(fc) = {ak)^e^^) follows the trend of 
the data remarkably well, while fitting functions with fewer free 
parameters (power law - f{k) = (ak)^) cannot model the shape 
accurately. 


More subtle effects such as the ‘beat-coupling’ (Hamilton 
et al. 2006; Rimes & Hamilton 2006; Sefusatti et al. 2006) 
and ‘local average’ (Sirko 2005; Takahashi et al. 2009; de 
Putter et al. 2012) will result in a cross-correlation even for 
large Ak. To account for those, we add a constant term to 
g{Ak). By combining the above effects we get 

g(Ak) = [ 7 V(Afc^ -b 7^)] + (1 - Q)sinc(u;Afc) -b /3 

The terms are combined in such a way as to enforce sr(0) = 1. 

Our final ansatz for the covariance matrix is 

Cij = PiPj f (ki) f {kj)g(ki — kj) (21) 

with functions f{k) and g{ki — kj) given by Eqs. (17) and 
(20) respectively. We will find the best-fitting numerical val¬ 
ues for free parameters a, b, v, a, 7 , uj, and /3, by fitting this 
ansatz to the sample covariance matrix measured from 600 
mocks. 

Figure 3 shows \/CV(/Pi (the fractional uncertainty in 
P{k) computed from the mocks) where the errorbars are 
computed using Eq. (10). Our fitting function provides an 
excellent fit to its shape. We find that the shape is difficult 
to recreate using functions with fewer free parameters, such 
as a simple power law. 

Figure 4 shows a similar plot for the off-diagonal ele¬ 
ments of the reduced covariance matrix, where we plotted 
the measurements in terms of Ak. Our fitting function seems 
to provide a good phenomenological description of this func¬ 
tion. We find that reducing the number of parameters by 
eliminating either the sine term (by setting q = 1 ) or the /3 
term (by setting /3 = 0 ) significantly worsens the fit. 

After performing the full fit to all independent covari¬ 
ance matrix elements we get a = (451 ± 35)h“^Mpc, 
b = -1.19 ± 0.02, u = (9.62 ± 0.32)h"^Mpc, 

a = 0.867 ±0.024, 7 = (5.17 ±0.16) x 10"^/i Mpc“\ 

u = (211.35 ± 6.14)/i“iMpc, and P = 0.0423 ± 0.0033, with 
= 250.6 for 269 degrees of freedom. The best-fitting 
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Figure 4. The independent ofE-diagonal elements of the correla¬ 
tion matrix plotted as a function of Afc. Grey data points repre¬ 
sent individual Vij elements while the black data points are the 
weighted mean and variance for the elements with the same value 
of ki — kj. Our proposed function, g{Ak), seems to fit quite well '1= 

for the full range of Afc (x^ = 219.6, Afjjop = 249). The best-fits <1 

for the Lorentzian plus constant term and Lorentzian plus sine 
term result in noticeably worse fits (x^ = 252.9, Afjjop = 251 and 
X^ = 388.7, Afpjop = 250, respectively). Note that the terms of 
each fitted function are combined in such a way as to enforce 
that they are equal to one when Afc = 0. 

value for ui is close to the theoretically expected value of 
the average depth of the survey divided by 27t. 


5 CONVERGENCE OF THE COVARIANCE 

MATRIX 

The main advantage of the fitting function approach is that 
the covariance matrix elements converge to their true val¬ 
ues much faster than the sample variance. Figure 5 shows 
the offset of individual covariance (and inverse covariance) 
matrix elements estimated using sample variance (light blue 
points) and our method (purple points) as the number of 
mocks increases. The offset is normalized to the standard 
deviation from the final {N = 600) covariance (and inverse 
covariance) matrix given by Eqs. (10) and (11) (or Eqs. (15) 
and (16) for the inverse). The fitting function method con¬ 
verges to the final result much faster both for the covariance 
and inverse covariance matrices. 

Figure 6 shows a histogram of the distribution of es¬ 
timated Cij elements around their true value for V = 50. 
This is a horizontal slice of the top panel of Figure 5. Al¬ 
ready, very few elements estimated with the fitting function 
method are outside Scr of the true covariance, while for the 
sample variance method the distribution is basically flat. 

At low N, there is a small bias in Cij elements deter¬ 
mined from the fitting function method but it’s significantly 
smaller than the variance in Cij-s determined from the sam¬ 
ple variance and therefore the fitting function method is still 
superior. 

Figure 7 shows the convergence to the final covariance 
in terms of per degree of freedom (DoF) for the two 
methods. Since the final covariance matrices computed with 
fitting function and sample variance methods are close to 



0 100 200 300 400 .300 600 


N (Number of Mocks) 

Figure 5. Convergence of the covariance (top) and inverse co- 
variance (bottom) matrices to their final values (computed using 
600 mocks) with the number of mocks. The light blue points show 
the matrices determined from the mock catalogs in the standard 
way, while the purple dots show the matrices as determined using 
our fitting function. It is clear that the matrices from the best¬ 
fitting functions converge much faster than the ones determined 
from sample variance. See the online article for the colour version 
of this plot. 

each other in x^, either one of them can be treated as a 
good approximation for the true covariance matrix. To put 
both methods on equal footing, when computing the 
for the fitting function method we compare it to the final 
{N = 600) sample variance method and vice versa. Figure 7 
clearly demonstrates that the fitting function method gen¬ 
erated Cij-s become consistent with the true covariance ma¬ 
trix much faster (already at A ~ 100) compared to those 
generated from the sample variance. 


6 CONCLUSIONS 

We propose a new method of estimating the elements of 
the covariance matrix for band-averaged P{k) measure¬ 
ments from galaxy surveys. The essence of the method is 
to find a fitting function for the covariance matrix and cal¬ 
ibrate its parameters using a sample covariance computed 
from a set of mock catalogues. We show that for the P{k) 
measurements from the BOSS DRll data in the range of 
0 ^ fc ^ 0.2h Mpc“^ a very simple, seven-parameter func- 
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Figure 6. Histograms of the differences between intermediate 
(generated using first 50 mocks) covariance matrices and the fi¬ 
nal (using all 600 mocks) covariance matrix normalized by the 
uncertainty in the final sample covariance matrix. While slightly 
biased, the fitting functions already provide a significantly better 
approximation of the final covariance matrix. 



0 100 200 300 400 500 600 

N (Number of Mocks) 


Figure 7. Reduced of tho covariance matrix compared to the 
final covariance. The solid line is for the results obtained using 
the fitting function and the dot-dashed line is for the sample co- 
variance matrix. Horizontal dashed lines show the expected 1, 2, 
and 3fT deviations from the x^/^DoF = 1 line. 


tion, given by Eq. (21), is sufficient to describe the covari¬ 
ance matrix. The fitting function covariance matrix is statis¬ 
tically indistinguishable from the sample covariance matrix 
computed from 600 mocks catalogues. 

The greatest advantage of this method, compared to 
the standard method of using the sample covariance ma¬ 
trix, is that it requires significantly fewer mock catalogues 
for calibration to converge to the true covariance matrix 
(see figure 5). For the BOSS DRll data, the fitting func¬ 
tion generated covariance matrices calibrated with as few 
as ~100 mocks were statistically indistinguishable from the 
sample covariance matrix generated with 600 mocks. To get 
a similar convergence with the sample covariance matrix we 
had to use more than 400 mocks (see figure 7). This advan¬ 
tage of the new method is especially relevant in situations 
where only a few mock catalogues are available. With only 


50 mock catalogues the distribution of the sample covari¬ 
ance around the true value is basically flat, while the fitting 
function method already provides a decent approximation 
(see figure 6). 

The specific functional form that we tested in this work 
may turn out to be an approximation that is too crude for 
future surveys, which will demand much higher accuracy on 
the determination of Cij. The functional form may have to 
be modified if one wishes to include P{k) measurements on 
much smaller scales as well. However, once an appropriate 
(for the sought precision) functional form is found, there 
is no reason why this method should not work with future 
data. 

In future work we plan to study how this new method 
works at higher precision. The uncertainties in the elements 
of the sample covariance (and inverse sample covariance) 
matrix scale inversely with the number of mocks (see se- 
ciont 3). A larger suite of mock catalogues would enable us 
to better estimate the true covariance (and inverse covari¬ 
ance) matrix with much smaller error bars on its elements. 
This would allow us to determine what modifications of our 
fitting function are required at higher precision to model the 
true covariance matrix. This would have the additional ben¬ 
efit of testing the fitting function method against a set of 
mocks which is completely independent from the set used in 
this work, ruling out the possibility that our successes were 
the result of some peculiarity which may be present in the 
data. 

This work was concerned only with the band-averaged 
P{k). Higher order Legendre polynomials of P{k) with re¬ 
spect to the line-of-sight provide valuable information about 
the nature of gravity and the expansion rate. The measure¬ 
ment of various Legendre moments of P{k) will be cross- 
correlated. We will address the question of modelling this 
larger covariance (and inverse covariance) matrix in future 
work as well. 
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